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ABSTRACT 


The  study  of  detonation  has  been  based  upon  hydrodynamic  theory. 
This  view  of  detonation  completely  ignores  the  actual  chemistry 
of  the  explosive  reaction.  Therefore,  the  dynamics  of  detonation 
on  a  molecular  level  remain  unknown.  The  purpose  of  this  project 
is  to  use  a  computer  to  investigate  the  propagation  of  detonation 
through  a  crystal. 

Research  in  this  area  is  hindered  by  the  fact  that  monitor¬ 
ing  instruments  are  destroyed  in  actual  detonations.  Computer 
simulations  avoid  this  problem  because  there  is  no  physical 
explosion.  The  actual  detonation  is  extremely  rapid;  collection 
of  data  at  designated  conditions  or  times  can  not  be  guaranteed. 
The  computer  does  not  have  this  problem  since  it  can  be  program¬ 
med  to  display  the  data  at  any  desired  condition  or  time. 

A  nonhomogeneous  crystal  of  diatomic  molecules  was  monitored 
to  discover  the  atomic  interactions  during  detonation.  A  Len- 
nard-Jones  potential  equation  was  used  to  represent  the  exo¬ 
thermic  reaction  between  diatomic  hydrogen  and  chlorine  mole¬ 
cules.  This  is  the  f i rst  project  to  use  the  natural  formation 
of  stable  reaction  products  to  achieve  exotherm i c 1 ty . 
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Computer  simulations  are  increasingly  used  in 
scientific  research.  They  can  be  used  to  verify  current 
theory  by  reproducing  laboratory  results  or  to  imitate  a 
system  which  can  not  be  empirically  monitored  in  the 
laboratory.  The  second  use  is  extremely  beneficial  because 
one  can  write  a  computer  program  which  generates  output  at 
any  given  time  or  condition.  The  computer  simulations  can 
isolate  particular  parameters  and  monitor  the  effects  of 
various  changes  to  these  parameters.  The  computer  simulation 
is  also  very  versatile,  since  it  can  focus  in  on  either 
macroscopic  or  microscopic  iews  of  a  system.  The  main 
advantage  of  a  computer  simulation  is  the  number  of 
computations  which  can  be  made  within  a  given  amount  of  time. 

My  project,  "A  Computer  Simulation  of  Detonation 
within  an  Energetic  Molecular  Crystal,”  took  advantage  of  the 
computer’s  ability  to  generate  and  monitor  a  system  which  can 
not  be  observed  through  ordinary  laboratory  means.  Problems 
with  monitoring  detonation  in  the  laboratory  are  that  any 
monitoring  devices  are  usually  destroyed  and  the  speed  of  the 
reaction  prevents  accurate  assimilation  of  data.  The 
computer  simulation  is  obviously  not  explosive  in  physical 
terms;  therefore,  the  first  problem  of  destroying  the 
monitoring  equipment  is  eliminated.  The  problem  of  the 
reaction’s  speed  is  easily  avoided  because  the  computer 
program  can  be  written  to  display  the  detonation  at  any  given 
time  or  condition  during  the  reaction.  The  data  generated  by 


the  computer  can  very  easily  be  stored  in  ordered  files  for 


future  use.  The  data  stored  is  also  selected  by  the 
programer.  He  can  generate  tables  of  numbers  or  even 
graphical  output.  These  advantages  make  computer  simulations 
a  very  good  medium  for  studying  detonation — if  the  proper 
equations  relating  the  interactions  between  the  atoms  can  be 
found.  The  main  focus  of  my  project  was  to  see  if  a  program 
could  solve  the  simultaneous  differential  equations  for  a 
large  number  of  particles  within  a  reasonable  amount  of 
computer  run  time. 

Previous  Uorks 

Scientists  have  done  studies  of  shock  waves  in  earlier 
simulations.  A.  M.  Karo,  J.  R.  Hardy,  and  F.  E.  Ualker  did 
work  with  a  homogeneous  monoatomic  crystal.  1  Their 
work  was  on  shock  induced  detonation.  They  used  Horse 
potentials  for  interaction  between  atoms,  but  the  atoms  only 
interacted  with  its  first  and  second  neighbors  in  the 
crystal.  If  the  crystal  array  is  disrupted,  the  atoms 
’’become  closely  juxtaposed  without  ’sensing’  one  another’s 
presence. "=  Although  this  greatly  reduces  the  run 
time  of  the  program,  it  allows  atoms  not  normally  next  to 
each  other  to  "pass  through"  one  another.  This  does  not  pose 
any  major  problems  when  low  values  of  initial  velocity  are 
given  to  the  crystal,  but  when  larger  values  of  velocity  or 
random  motion  perpendicular  to  the  shock  wave  propagation  are 
used  the  crystal  structure  disintegrates.3 
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Karo,  Hardy,  and  Ualker  also  observed  that  part  of  the 
crystal  departed  the  main  lattice  structure  when  a  detonation 
wave  impinged  a  free  surface."*  These  problems  can  be 
solved  if  they  allow  each  atom  to  interact  with  every  other 
atom  in  the  crystal  and  if  they  reduce  the  number  of  atoms 
located  on  a  free  surface  in  the  crystal  being  observed. 

Two  types  of  potential  equations  were  used:  endo¬ 
thermic  and  exothermic.  These  equations  were  used  during 
separate  computer  runs.  The  endothermic  potential  equation 
represented  the  breaking  of  bonds  between  adjacent  atoms. 

Th  runs  made  with  this  potential  function  had  a  quiescent 
behavior.  The  exothermic  potential  equation  was  written  by 
the  authors  to  be  an  exothermic  response  to  a  bond  stretched 
beyond  a  certain  limit.  This  reaction  does  not  occur  in 
nature.  The  entire  lattice  was  torn  apart  when  this  unnatural 
exothermic  equation  was  used.®  These  equations  are 
fine  for  initial  research  but  they  are  not  the  results  of  a 
natural  detonation.  The  energy  from  detonation  comes  from 
atoms,  molecules,  and  molecular  fragments  forming  more  stable 
molecules  after  the  detonation  wave  has  passed.  The  energy 
released  by  the  more  stable  molecules  is  what  propagates  the 
detonation  wa>'e  and  causes  the  thermal  and  shock  effects  of 
the  explosion. 

D.  H  Tsai  and  S.  F.  Trevino  did  studies  on  shock  wave 
propagation  through  a  homogeneous  diatomic  crystal.  This 
work  is  closer  to  actual  detonation,  because  their 


hypothetical  diatomic  molecules  can  break  apart  and  remain 


unattached  or  they  can  reform  diatomic  molecules.  The 
reformed  molecule  is  the  same  molecule  which  was  originally 
present;  therefore  more  stable  substances  are  not  produced. 

A  potential  energy  equation  which  is  nonoccurring  in  nature 
must  be  used  to  make  the  equation  exothermic.  They  used  two 
different  Horse  potential  equations:  one  for  dissociated 

atoms  and  one  for  atoms  bonded  together  as  a  diatomic 
molecule.  The  dissociation  of  the  atoms  was  given  an 
exotherm  1 c i ty  factor,  like  Karo,  Hardy,  and  Ualker  used, 
which  does  not  occur  in  nature.  Each  atom  had  a  flag  to  show 
whether  or  not  it  was  bonded  to  another  atom.  The  value  of 
this  flag  would  determine  which  of  the  two  potential  energy 
equations  would  be  used.  Figures  1  and  2  are  of  the 
i nt r amo 1  ecu  1 ar  and  l ntermo 1  ecu  1 ar  potentials  respectively. 

Note  the  strong  i ntermo 1  ecu  1 ar  potential  in  Figure  2;  normally 
i ntr amo 1  ecu  1 ar  forces  are  stronger  than  i ntermo 1  ecu  1 ar  ones. 
Figure  3  shows  the  endothermic  and  exothermic  potential  curves 
used  by  Tsai  and  Trevino.  The  upper  curve  is  the 
intramolecular  potential;  it  does  not  take  much  activation 
energy  to  get  it  to  reach  the  l ntermo 1  ecu  1 ar  curve.  Once  an 
atom  vs  "influenced”  by  the  i ntermo 1  ecu  1 ar  curve,  it  will 
dissociate  to  become  stable.  The  stability  gives  the  reaction 
energy,  but  as  mentioned  before,  it  is  not  a  natural 
phenomenon  .  & 


R.  A.  riacDonald  and  D.  H.  Tsai  worked  on  dynamical 
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Figure  1 

Intramolecular  potential  equation  invented  by  Tsai  and 

Trevino 


Energy  units  are  arbirary. 


Distance  bstwssn  atoms 


Figure  2 

Intermolecular  potential  equation  invented  by  Tsai  and  Trevino 
Energy  units  are  arbitrary. 


Distance  between  atoms 


Figure  3 

The  change  in  potential  from  intramolecular  to  intermolecular 
interactions.  The  energy  units  are  arbitrary. 
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calculations  of  energy  transport  in  crystalline  solids. 

Their  study  uas  testing  whether  or  not  thermal  equilibrium 
existed  behind  the  shock  front  in  a  solid.  They  concluded 
the  thermal  equilibrium  propagated  at  such  a  slow  velocity 
that  thermal  equilibrium  uas  not  achieved  behind  the  shock 
wave.  The  computed  results  showed  that  a  relaxation  region 
existed  behind  the  shock  wave,  but  that  it  never  reached 
equilibrium.  The  kinetic  and  potential  energies  would  not 
transmit  the  equilibrium  fast  enough  to  follow  the  shock 
wave.  Their  results  did  show  that  the  relaxation  time  is 
dependent  upon  the  internal  degrees  of  freedom  of  the 
molecule.  They  originally  only  had  included  coupling  between 
molecules  and  not  i ntramo 1  ecu  1 ar  degrees  of  freedom.  Uhen 
the  intramolecular  degrees  of  freedom  were  included,  the 
relaxation  time  was  shortened.7. 

Initial  Conditions 

To  avoid  some  of  the  problems  encountered  by  the  other 
scientists,  I  utilized  several  different  initial  conditions. 
These  differences  included  original  atomic  spacing,  periodic 
boundaries,  expansion  of  the  crystal,  and  interaction  between 
all  atoms.  The  crystal  was  not  set  up  starting  at  the 
origin.  Instead  the  crystal  was  shifted  out  along  the  x-axis 
15  Angstroms  and  shifted  up  one  half  of  the  equilibrium 
spacing  between  molecules  along  the  y-axis.  (Figure  4)  The 
reason  for  moving  out  in  the  x  direction  is  to  allow  the 
atoms  receiving  initial  momentum  and  velocity  to  approach  the 
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Monoatomic  crystal  generated  by  computer  in  perfect  rows  and 
columns.  Note  that  monoatomic  molecules  are  used  for  simplicity 
of  showing  the  initial  parameters. 


Figure  4 
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crystal  with  natural  interactions.  It  is  not  correct  for 
matter  to  instantaneously  change  states.  (Even  capacitors 
take  a  small  amount  of  time  to  charge  up  in  electronic 
circuitry).  The  rest  of  the  crystal  is  also  allowed  to 
settle  into  equilibrium  position  before  the  detonation  is 
initiated.  The  shift  in  the  y  direction  is  related  to  the 
periodic  boundaries  and  the  expansion  of  the  crystal. 

Periodic  boundaries  will  eliminate  the  problem  of 
hawing  large  numbers  of  observed  atoms  along  free  surfaces. 
The  top  and  bottom  edges  make  up  the  majority  of  the 
crystal’s  free  surface  area.  Approx i mate  1 y  forty  percent  of 
all  atoms  involved  in  the  simulation  are  situated  along  the 
edge.  This  would  be  an  unrealistic  ratio  in  an  actual 
physical  crystal.  Physical  crystals  have  most  of  their  atoms 
on  the  inside,  while  only  a  small  ratio  of  the  atoms  are 
located  on  a  free  surface.  There  are  two  ways  to  solve  this 
problem:  use  several  million  atoms  in  the  simulation  or 

expand  the  crystal.  The  use  of  millions  of  atoms  would  make 
the  program’s  run  time  on  the  order  of  years  instead  of  days 
at  the  current  state  of  computer  technology.  Using  that 
large  a  number  of  atoms  in  a  simulation  is  not  feasible.  The 
use  of  periodic  boundaries  only  doubles  the  number  of  atoms 
involved  in  the  simulation,  and  significantly  lowers  the 
ratio  of  atoms  along  a  free  surface  to  approximately  four 
percent.  Periodic  boundaries  were  set  up  to  keep  the  atoms 


contained  in  the  crystalline  area 


This  avoids  the  problem 


experienced  by  Karo,  Hardy,  and  Ualker  in  which  the  atoms 
along  a  free  surface  departed  the  crystal  structure  into  free 
space.  fly  periodic  boundaries  uere  set  at  one  half  the 
equilibrium  distance  between  molecules  in  the  y  direction. 
This  shifts  the  crystal  up  one  half  the  equilibrium  distance 
along  the  y-axis.  (Figure  5)  The  other  periodic  boundary 
exists  one  half  the  equilibrium  distance  above  the  center  of 
the  initial  position  of  the  top  row  of  atoms.  These 
boundaries  remain  in  the  same  position  throughout  the 
computer  run.  If  an  atom  ever  crosses  a  boundary,  (Figure 
6),  then  it  is  shifted  to  the  other  periodic  boundary. 

(Figure  7)  The  shifted  atom  only  changes  position  in  the 
y-coord i nate ;  the  x-coord i nate ,  the  velocity,  the  momentum, 
and  the  momentum  derivative  remain  the  same.  If  the  atom  in 
Figure  6  has  a  y-coor i d i nate  .002  Angstroms  below  the 
periodic  boundary  at  the  y-origin,  then  its  y-coordinate  is 
shifted  to  .002  Angstroms  below  the  periodic  boundary  on  the 
top.  This  effectively  wraps  the  crystal  around  along  the 
y-coor l d l nates .  (Figure  8)  The  expansion  of  the  crystal  is 
an  extension  of  the  periodic  boundaries;  the  expansion 
enables  the  movement  of  atoms  from  one  periodic  boundary  to 
the  other  to  be  a  natural  occurrence. 

Expanding  the  crystal  involves  generating  a  new  set  of 
atoms  by  taking  atoms  with  y-coord i nates  situated  between 
the  origin  and  one  half  of  the  top  periodic  boundary,  (Figure 
9),  and  shifting  them  up  above  the  top  periodic  boundary. 


„  BflP-odic  boundary 
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Positions  of  the  periodic  boundaries  around  the  crystal 


Figure  15 
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Atom  crossing  lower  periodic  boundary 


Figure  6 
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Atom  shifted  to  top  of  crystal  after  it  has  crossed  periodic  boundary 


Figure  7 
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Atoms  with  shading  will  be  projected  abouve  the  crystal  during 
"crystal  expansion." 


Figure  9 
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(Figure  10)  Then  the  remaining  atoms,  those  in  the  top  half 
of  the  original  crystal,  are  shifted  belou  the  bottom 
periodic  boundary.  This  has  doubled  the  number  of  atoms 
involved  in  the  simulation.  (Figure  11)  The  neu  atoms  are 
contained  only  in  temporary  memory,  and  they  are  regenerated 
during  each  time  step.  They  only  have  positions  both  a  and 
y;  they  do  not  have  their  oun  velocity,  momentum,  or  momentum 
derivative.  These  parameters  are  not  needed  since  the  neu 
atoms  are  projected  images  of  the  original  atoms.  The 
velocity,  momentum,  and  momentum  derivatives  from  the 
original  atoms  uill  correctly  move  the  neu  expanded  atoms. 

The  reason  for  having  the  neu  expanded  atoms  is  to  enclose 
the  original  crystal.  The  atoms  from  the  original  crystal 
are  nou  interior  atoms,  uhich  uill  give  more  accurate 
predictions  than  atoms  impinging  a  free  surface  edge. 

Looking  back  at  Figure  6,  there  is  a  ’’hole”  generated  above 
the  column  uhich  is  shifted  doun .  By  expanding  the  crystal 
this  hole  has  been  filled  by  an  expanded  atom  since  the 
beginning  of  the  program.  Nou  that  the  original  atom  has 
crossed  the  periodic  boundary,  the  original  atom  is  moved 
into  this  hole.  One  might  think  that  this  movement  creates  a 
neu  hole  uhere  the  original  atom  uas  moved  from.  This  is  not 
the  case  because  the  expanded  crystal  nou  generates  an  atom 
in  this  hole  during  subsequent  time  steps.  The  expanded 
atoms,  therefore,  are  actively  interacting  uith  the  original 
atoms  during  each  time  step. 
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Projected  atoms  generated  during  crystal  expansion 


Figure  10 
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Crystal  generated  after  expansion.  Note  that 
the  space  in  between  the  new  atoms  and  the  or¬ 
iginal  crystal  was  artificially  put  in  the  figure 
so  that  the  difference  is  easily  noticed. 


Figure  II 
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Potential  Equations 


The  interaction  between  the  atoms  is  caused  by  the 
potential  equations.  Two  potential  equations  are  used  in 
computer  simulations:  the  Horse  equation  and  the 

Lennard-Jones  equation.  The  two  are  very  similar — except  for 
the  behavior  of  the  atoms  at  small  interatomic  distances. 

The  Horse  equation  has  a  finite  potential  value  when  the 
interatomic  distance  goes  to  zero.  This  accounts  for  atoms 
"passing  through”  each  other  which  is  physically  impossible. 
Karo,  Hardy,  and  Ualker  had  stated  that  the  phenomenon  of 
atoms  passing  through  one  another  occurred  during  their 
simulations.  They  had  used  the  Horse  potential  equation 
which  contributed  to  this  error.  The  Lennard-Jones  potential 
equation  goes  to  infinity  as  the  interatomic  distance  goes  to 
zero.  (Figures  12-14)  Atoms  experience  this  repulsion  in  real 
life;  therefore,  the  Lennard=Jones  equation  is  more  accurate 
when  close  interatomic  distances  will  be  experienced.  The 
Lennard-Jones  equation  is  normally  given  in  the  form: 

_i  r  i  s  4  *  *  •  O  I"  1  1  -  -•  1  O  r  1 

r  -  inter-atomic  distance 

Lennard-Jones  potential  equations  are  normally  inter- 
molecular.  For  this  simulation  the  parameters  were  changed 
to  have  intramolecular  interactions.  This  was  done  to  allow 


an  atom  to  interact  within  the  original  molecule  it  was 
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Interatomic  Distance  (Angstroms) 


Figure  13 

Lennard-Jones  potential  for  interatomic  H-Cl  interactions 
The  energy  units  are  amu  R^/  psec^. 
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assigned  to,  while  being  able  to  interact  uith  the 
surrounding  atoms  to  form  the  most  stable  compound  possible. 
This  uas  simply  done  by  taking  the  equilibrium  1 ntramo 1  ecu  1 ar 


distance  for  each  of  the  three  molecules,  HH ,  HC1 ,  and  C1C1, 


and  solving  for  the  constants  o  and 


The  equ i 1 i br i urn 


interatomic  distance  is  related  to  the  o  term  by  a  multiple 
of  21-'*.  The  term  is  a  measure  of  the  extent  of 
attraction  between  pairs  of  molecules.  The  values  for  the 
attraction  uas  taken  from  flolecular  Theory  of  Gases  and 
Liquids  by  J.  0.  H i rschf e 1 der ,  C.  F.  Curtiss,  and  R.  B. 

Bird.13  Units  were  converted  to  a  system  of 

Angstroms,  picoseconds,  and  atomic  mass  units  to  avoid  the 
use  of  large  exponents  uith  the  potential  energy  values. 

This  system  also  made  the  position  values  manageable  numbers 
uith  the  Angstrom  unit  measure.  Lennard- Jones  potential 
equations  are  very  accurate  for  non-polar  molecules,  HH  and 
C1C1;  however  for  polar  molecules,  HC1 ,  the  equations  are  not 
as  accurate.  The  Lennard- Jones  potential  "may  be  useful  for 
purposes  of  calculations  until  the  theory  needed  for 
describing  complex  molecules  has  been  developed.”'” 

The  Lennard- Jones  potential  equations;  therefore,  are 
best  for  the  requirements  of  my  simulation.  To  save 
calculation  time  the  form  of  the  equation  uas  converted  to: 


8  M  (  1  /  r  i  1 : 


*  <  t  r  > ■ 


The  constants  4,  ,  and  o  were  raised  to  the  appropriate 

power  and  then  multiplied  together.  This  alteration  saves 


six  calculations  per  atom-atom  interaction.  The  interatomic 
distance  raised  to  the  tuelth  pouer  causes  the  repulsion 
betueen  atoms  at  close  distances.  (Figure  12)  The  negative  of 
the  sixth  pouer  of  the  interatomic  distance  produces  the 
bottom  of  the  well,  so  that  there  is  an  equilibrium  distance 
betueen  the  atoms.10 

Equations  of  Notion 

The  equations  of  motion  are  the  momenta  and  coordinate 
equations  for  each  atom.  These  exist  in  all  three 
coordinates,  but  my  computer  simulation  is  limited  to  tuo 
dimensions.  The  Hamiltonian  equations  are: 

q*  =  JH  /  d pi. 

=  J  T  /  d  p  * 
p*  =  -  JH  /  d  q»- 

—  -•'i'd  >  <J  q  u 

During  each  time  step  of  the  simulation,  each  atom  has  a 
set  of  these  equations  uith  every  other  atom  in  the  expanded 
crystal.  The  momenta  and  coordinate  equations  are  directly 
dependent  upon  one  another;  they  are  simultaneous 
differential  equations.  The  kinetic  energy,  T,  is  equal  to 
the  summation  of  the  momentum  squared  divided  by  the  tuice 
the  mass  of  the  atom.  The  V  term,  the  potential  energy,  is 
the  Lennard- Jones  equation.  The  identity: 

J-F  '  d q*  -  ‘di/dr)  idr/dq*> 

is  used  to  solve  the  momenta  equations,  since  the  Lennard- 
Jones  equation  must  be  solved  uith  respect  to  the  coordinate 
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derivative.  The  equations  are: 

T  —  p  *  (  2  *  iTl  fe  < 

C|1  k  —  pi  U  /  1  flt-r 

V  =  B*r-i:2-A*r 
r  j  k  ~  (  q  i  j  ~  q  1 1< )  *  1  ' 

pin  —  -dV  /  Jqik 

—  (  —  o'V  /  dr  >  (dr  /  d q i  _i  > 

subscript  i  refers  to  x  or  y-coordinate 
These  equations  are  in  a  usable  form  for  solving  the 
atomic  interactions.  Substituting  in  the  actual 
Lennard-Jones  equation  for  pin  gives  the  equation  used 
in  the  subroutine  which  calculates  the  force  on  each  atom. 
This  equation  is: 

pm  =  ( 1 1  *  +  o*A*r~°>  * 

(qik  -  q  an  i  - 1  i  1 1 

These  equations  are  now  ready  to  be  integrated  to  calculate 
the  momenta  and  coordinates  for  each  atom. 

Integration  Techniques 

Now  let  us  address  the  problem  of  properly  integrating 
the  equations  of  motion  and  momentum.  The  main  problem  is 
that  the  program  has  several  simultaneous  differential 
equations  for  each  atom.  This  requires  five  hundred  seventy- 
six  sets  of  simultaneous  equations  to  be  solved  during  each 
time  step.  Ilany  pre-ur  itten  integration  programs,  such  as 
DGEAR  from  the  International  Mathematics  and  Science  Library, 
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ItlSL,  can  not  efficiently  handle  such  a  large  number  of 
simultaneous  equations.  Research  into  integration  techniques 
determined  that  the  Euler  methods  might  be  uork. 

Unfortunately  the  Euler  method  generates  an  error  in  the 
calculation  of  a  new  position,  Xn»i,  if  the  velocity 
is  not  constant  over  the  time  interval.  The  Euler  method  is 
a  most  simple  method,  but  it  lacks  accuracy.  An  extremely 
small  step  size  is  required  to  get  any  accuracy.  11 
In  real  life  the  atoms  interact  at  infinitely  small  time 
intervals,  i.e.  continuous  integration.  Computer  simulations 
require  that  a  reasonable  time  interval  be  selected;  since 
the  forces  acting  upon  each  atom  changes  with  its  position, 
the  velocity  is  constantly  changing.  This  generates  an  error 
over  time  increments  in  which  the  velocity  is  not  constant. 
Each  integration  with  the  Euler  method  has  an  error.  The 
error  starts  out  small,  but  each  subsequent  integration 
increases  the  error  from  the  correct  value.  The  Euler  method 
becomes  more  accurate  the  smaller  the  time  increment,  but  the 
error  generated  by  a  time  increment  which  has  a  reasonable 
program  run  time  is  unacceptably  large.  Fortunately  there  is 
a  Modified  Euler  method  which  is  over  twenty-two  times  more 
accurate  for  each  integration  over  the  same  time 
increment.  1  ,;2  It  should  be  noted  that  as  the  time 
increment  becomes  infinitely  small,  both  the  Euler  method  and 
the  Modified  Euler  method  will  give  the  same  correct  answer. 
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The  Modified  Euler  method  uses  the  arithmetic  average  of  the 
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velocities  at  the  beginning  and  end  of  the  time  increment. 
This  average  velocity  more  accurately  predicts  the  true 
position.  The  equation  for  the  Modified  Euler  method  is: 

X r%  —  1  —  Xn  .  5  ^  (  Vn  r  'Xr,  -*■  X  1  #  tl 

X  -  position  V  -  velocity  h  -  time  increment 
One  can  not  use  the  Modified  Euler  method  directly  from  the 
equation.  The  equation  requires  that  Un»i  be  known 
to  predict  the  value  of  Xn«i.  The  Euler  method  is 
used  to  predict  a  value  for  Xo*t;  this  value  is 
temporarily  stored  to  calculate  a  value  for  Un»i. 

This  value  of  Vn*i  is  then  used  in  the  Modified 
Euler  equation  to  predict  an  accurate  value  for  Xr->~i. 

Other  integration  methods  exist,  example  the  Runge-Kutta 
method,  but  their  run  time  makes  them  inefficient  for  this 
particular  problem.  The  Runge-Kutta  method  is  a  fourth  order 
solution,  which  doubles  the  run  time  and  storage  required  for 
the  program. 

The  Program 

The  actual  program  is  included  in  Appendix  A.  It  is 
written  in  Fortran  77.  The  program  starts  out  by  dimension¬ 
ing  all  of  the  arrays  required  for  the  program.  The 
position,  velocity,  momentum,  and  momentum  derivative  are  set 
up  for  the  x  and  y  terms  for  each  atom  in  the  crystal.  The  x 
and  y  position  for  the  expanded  crystal  are  also  included.  A 
second  set  of  arrays  for  position,  velocity,  momentum,  and 
momentum  derivative  has  been  set  up  to  use  the  Modified  Euler 


integration  method.  The  final  array  dimensioned  is  to  record 
the  type  (T)  of  the  atom.  The  parameters  for  the  crystal 
size  follow  the  arrays.  The  LENGTH  is  a  summation  of  the 
equilibrium  distance  between  two  hydrogen  atoms,  a  HH 
molecule  and  a  C1C1  molecule,  two  chlorine  atoms,  and  a  HH 
molecule  and  a  C1C1  molecule  again.  This  is  the  periodic 
cycle  to  set  up  the  initial  crystal.  The  equilibrium 
distance  between  a  HH  molecule  and  a  C1C1  molecule  was 
received  from  A.  Blumen  and  C.  Merkel  in  their  ab 
initio  study  of  hydrogen  chloride.  13  The 
YI10VE  term  is  the  top  periodic  boundary. 

The  opening  of  external  files  is  the  next  major 
portion  of  the  program.  The  ’ totena ’  file  stores  the  total 
energy  of  the  crystal;  the  total  energy  of  the  crystal  is 
periodically  calculated  to  monitor  the  validity  of  the 
simulation.  If  the  simulation  is  bad,  then  the  total  energy 
will  exponentially  increase.  The  'savfla'  file  saves  the 
type,  position,  and  momentum  of  each  atom  at  designated 
times,  so  that  if  the  computer  crashes  the  program  can  be 
continued  without  restarting  from  the  beginning.  The  ’  ta  ’ 
and  ’ ia 1  files  store  stop  action  frames  of  the  crystal. 

These  files  are  generated  to  display  the  type,  size,  and 
position  of  each  of  the  original  atoms  in  the  crystal. 

The  subroutines  which  set  up  the  initial  crystal  and  give 
initial  momentum  to  selected  atoms  are  called.  These 


subroutines  are  towards  the  end  of  the  program.  The  POSition 
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subroutine  starts  on  page  A-9.  The  initial  positions  are 
constructed  by  the  use  of  two  Do  loops;  one  is  used  for  rows 
and  the  other  is  for  columns.  As  -ach  atom  is  assigned  a 
position,  it  is  also  assigned  a  Type.  The  rows  alternate 
uhich  type  of  atom  is  on  the  left  side  of  the  crystal.  The 
alternation  is  flagged  by  the  value  of  the  ’J’  variable.  To 
vertically  align  the  center  of  each  of  the  diatomic  molecules 
the  ’DISTCENT’  term  is  used  on  every  other  row.  Finally  the 
first  four  atoms  in  each  rou  are  moved  back  to  allow  the 
crystal  to  interact  naturally. 

The  norientum  subroutine  is  located  on  page  A-?.  This 
subroutine  sets  all  initial  momenta  to  zero;  then  it  gives 
the  first  four  atoms  in  each  rou  initial  momentum.  The 
hydrogen  and  chlorine  atoms  get  different  initial  momentum 
values.  The  values  of  the  momentum  in  the  x  direction  are 
calculated  so  that  the  atoms  get  the  same  initial  velocity. 
The  momentum  in  the  y  direction  is  randomly  generated  to 
simulate  random  thermal  energy. 

A  comment  block  stating  the  units,  back  in  the  main 
program,  follows  the  return  from  the  non  subroutine.  The 
initial  time,  time  increment,  and  final  time  are  now  set. 
These  three  terms  are  used  in  a  Do  loop,  which  contains  the 
integration  and  prediction  portions  of  the  Hod  1 f i ed  Euler 
method.  The  initial  positions  of  the  atoms  are  sent  to  a 
file  through  the  URITFILE  subroutine.  The  URITFILE 
subroutine  is  on  page  A-ll.  This  subroutine  generates  a  file 
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which  displays  the  positions  of  the  atoms  in  the  crystal.  An 
abridged  output  file  is  included  in  Appendix  A.  The  middle 
atoms  were  omitted  from  the  appendix  to  keep  the  data  at  a 
reasonable  level.  The  data  file  contains  commands  for  the 
PS300  computer  to  display  and  annotate  the  crystal.  Each 
atom  must  be  stored  as  a  list  of  vectors,  so  that  it  can  be 
displayed  in  a  circular  form.  The  atom  type  determines  which 
color  the  display  is.  Hydrogen  atoms  are  red  and  chlorine 
are  green;  these  are  the  standard  colors  used  to  represent 
these  two  atoms.  The  data  stored  by  the  URITFILE  subroutine 
will  be  used  to  animate  the  molecular  interactions. 

The  LJ  subroutine  is  the  next  major  portion  of  the 
program.  LJ  stands  for  Lennard- Jones ,  because  this 
subroutine  calculates  the  Lennard- Jones  potential  between 
atoms.  The  LJ  subroutine  is  listed  on  page  A-5 .  The 
momentum  derivatives  are  reset  to  zero,  so  that  the  summation 
of  each  atom’s  total  potential  energy  is  only  for  the  current 
time  step.  The  velocity  is  then  calculated  by  dividing  each 
atom’s  momentum  by  the  appropriate  mass.  Expansion  of  the 
crystal  is  the  next  major  routine.  This  is  done  by  comparing 
the  y-coordinate  position  with  the  center  of  the  crystal.  If 
the  atom  is  above  the  center  of  the  crystal,  then  another 
atom  is  generated  beneath  the  bottom  periodic  boundary  the 
same  distance  the  original  atom  was  beneath  the  top  periodic 
boundary.  Both  of  these  atoms  have  the  same  x-coordinate 
value.  If  the  original  atom  is  below  the  center  of  the 
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crystal,  then  an  expanded  crystal  atom  vs  generated  above  the 
top  periodic  boundary.  Now  the  subroutine  AL  is  called  to 
calculate  the  potential  forces  between  the  atoms. 

Subroutine  AL  simply  calculates  the  potential  force 
for  each  atom  in  the  original  crystal.  The  distance  between 
the  atoms  in  both  the  x  and  y-coord i nates  vs  calculated,  then 
these  two  terms  are  combined  to  find  the  magnitude  of  the 
vector  distance  between  the  two  atoms.  This  distance  is  then 
used  in  the  derivative  of  the  Lennard- Jones  potential  energy 
equation  to  calculate  the  force  between  the  atoms.  The  force 
is  stored  in  the  momentum  derivative  array.  The  computer 
determines  whether  the  atoms  are  two  hydrogens,  two 
chlorines,  or  a  combination  so  that  the  proper  coefficients 
may  be  used.  The  direction  of  the  force  is  determined  by  the 
magnitude,  positive  or  negative,  of  the  x  and  y  distances. 

The  subroutine  moves  on  to  the  next  atom  of  the  original 
crystal  when  it  has  calculated  the  force  between  this  atom 
and  all  other  atoms  within  the  expanded  crystal.  It  has  been 
determined  that  interatomic  distances  of  greater  than  ten 
Angstroms  will  not  generate  a  significant  force.  Therefore, 
if  the  magnitude  of  the  vector  distance  is  greater  than  ten, 
then  the  calculation  of  the  force  is  skipped  to  save  run 
time.  After  all  of  the  original  atoms  have  had  their  new 
forces  calculated  the  AL  subroutine  returns  to  the  LJ 
subroutine.  The  LJ  subroutine  immediately  returns  to  the  main 
program  to  be  used  in  the  Modified  Euler  integration  method. 
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The  velocity  and  momentum  derivative  values  returned 
from  the  LJ  subroutine  are  used  to  predict  the  position  and 
momentum  values  by  the  Euler  method.  These  new  values  are 
sent  back  to  the  LJ  subroutine  to  calculate  a  neu  set  of 
velocities  and  momentum  derivatives.  The  LJ  subroutine 
repeats  the  force  calculation  process,  and  returns  the  second 
set  of  velocities  and  momentum  derivatives.  The  first  and 
second  sets  of  velocities  and  momentum  drivatives  are  used  in 
the  Ilodified  Euler  integration  method  to  calculate  the 
correct  neu  positions  and  momenta. 

The  program  nou  checks  to  see  if  the  current  time  is 
one  of  the  designated  output  flags.  If  it  is  then  an  output 
file  is  generated,  the  total  energy  of  the  system  is 
calculated,  and  the  important  parameter  data  is  stored. 
Finally  the  original  atoms  are  checked  to  see  if  they  have 
crossed  a  periodic  boundary.  If  one  has,  it  is  moved  to  its 
appro;.  Mate  position  near  the  other  periodic  boundary.  This 
completes  one  pass  for  the  computer  run;  the  program  will  nou 
loop  through  this  routine  until  the  final  time  is  reached . 

Cone  1  us  1  ons 

The  program  is  the  first  computer  simulation  of 
detonation  using  either  a  realistic  exothermic  reaction  or 
nonhomogeneous  diatomic  molecules.  This  is  a  major  step 
touards  simulations  of  explosives  which  are  in  current 
arsenals.  These  molecules  are  too  large  to  have  a  computer 


simulation  with  present  technology. 


I  used  the  Lennard- Jones 


potential  equation  for  a  hydrogen  and  chlorine  crystal. 
(Figure  15)  The  original  crystal  consisted  of  HH  and  C1C1 
molecules,  as  the  detonation  progressed  the  molecules 
separated  and  began  to  form  HC 1 .  (Figure  16) 

The  run  time  on  the  DEC  UAX  11/780  computer  is 
approx i mate  1 y  six  days  to  simulate  one  picosecond.  This  is 
actual  working  CPU  time — not  elapsed  clock  time.'  The  speed 
of  the  UAX  is  one  million  instructions  per  second,  1  mips.  A 
supercomputer  with  the  speed  of  twenty  mips  could  easily 
handle  this  program.  Because  of  the  length  of  run  time  and 
the  periodic  computer  crashes,  the  amount  of  data  I  could 
obtain  is  limited. 

The  data  for  each  simulation  run  is  collected  in  data 
files;  an  example  is  shown  on  pages  A-14  to  A-20 .  These 
files  display  a  color  representat i on  of  the  crystal.  The 
crystal  is  two  dimensional  and  the  display  accurately  shows 
the  progression  of  the  detonation  wave.  The  only  drawback  to 
my  program  is  that  it  allows  the  atoms  to  clump  together. 
Since  this  is  only  a  simulation,  the  clumping  factor  is  not 
very  detrimental.  In  fact  it  increases  the  exotherm l c l ty  of 
the  react i on . 

Future  work  made  possible  by  my  project  will  explore 
different  parameters  of  detonation.  Changing  the  crystal 
orientation,  different  exotherm i c i ty  values  of  reaction, 
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A- la 


c** 

★  * 

c** 

★  ★ 

c** 

LENNARD- JONES  (6-12)  POTENTIAL 

*  * 

C** 

WITH  HH  &  C1C1  MOLECULES  IN  TWO  DIMENSIONS 

★  ★ 

c** 

★  ★ 

c** 

by 

★  * 

c** 

Alan  D.  Boyd 

** 

c** 

Midshipman  First  Class 

' 

*  * 

c** 

★  * 

c** 

Trident  Project  '86 

*  * 

c** 

★  ★ 

DIMENSION  ALL  ARRAYS  REQUIRED 


DIMENSION  X (288) , Y (288) 

DIMENSION  XM0M(144) ,  YM0M(144) 
DIMENSION  XVEL  (144)  ,  YVEL  (144) 
DIMENSION  XMOMD  (144)  ,  YMOMD  (144) 
DIMENSION  X2  (288)  ,Y2(288) 
DIMENSION  XMOM2 (144) , YM0M2 (144) 
DIMENSION  XVEL2 (144) ,YVEL2(144) 
DIMENSION  XMOMD2 (144) , YMOMD2 (144) 
DIMENSION  T (144) 


*f  *  ★  ★ 

*t  *  *  ★ 

**  SET  UP  ALL  COMMON  VARIABLES  ** 

**  ** 

*  *  *  * 


COMMON/B1/ORIGMOM 

*  it  *  * 

**  SET  UP  PARAMETERS  FOR  CRYSTAL  SIZE  ** 

★  ★  ★  * 


PARAMETER 
Check  value 
PARAMETER 
PARAMETER 
PARAMETER 


(XORIGIN  =15.) 

for  sizes  and  equation  for  ymove 
(LENGTH  =  10.11411) 

(SIZES  =  3.69) 

(YMOVE  =  6 . *SIZES) 


THE  NUMBER  OF  ATOMS  IN  THE  (ORIGINAL)  CRYSTAL  IS  =  N 


v.-z.v 


.N'.'.’-.'VvViV.V 


n  n  no 
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PARAMETER  (N  =  144) 
C 
C 


c** 

★  * 

c** 

OPEN  FILES  FOR  OUTPUT 

★  ★ 

c** 

** 

c** 

FIRST  ONE  IS  FOR  CALCULATION 

OF  TOTAL  ENERGY 

*  * 

c** 

SECOND  ONE  IS  FOR  STORAGE  IN 

CASE  OF  SHUTDOWN 

★  * 

c** 

THE  REST  ARE  FOR  STORING  THE 

ATOMS'  POSITIONS 

** 

c** 

AT  VARIOUS  INCREMENTS  OF  TFINAL 

★  * 

c** 

★  * 

OPEN  (UNIT=1 ,  STATUS="  UNKNOWN" ,  FILE="totena") 
OPEN (UNIT=2 , STATUS="UNKNOWN" , FILE="savf la") 
OPEN (UNIT=3  , STATUS= "UNKNOWN" , FILE="taO . 300") 
OPEN (UNIT=4  , STATUS=" UNKNOWN" , FILE="ial . 300") 
OPEN (UNIT=10, STATUS="UNKNOWN" , FILE="tal .300") 
OPEN (UNIT=11, STATUS="UNKNOWN",FILE="ta2 . 300") 
OPEN (UNIT=12 , STATUS=" UNKNOWN" , EILE="ta3 .300") 
OPEN (UNIT=13, STATUS="UNKNOWN" , FILE="ta4 . 300") 
OPEN (UNIT=14, STATUS="UNKNOWN", FILE="ta5 . 300") 
OPEN  (UNIT=15/  STATUS="UNKNOWN" ,  FILE="t:a6 .300") 
OPEN (UNIT=16/ STATUS="UNKNOWN" , FILE="ta7 . 300") 
OPEN (UNIT=17, STATUS=" UNKNOWN" , FILE="ta8 . 300") 
OPEN (UNIT=18/ STATUS="UNKN0WN".FILE="ta9 .300") 
OPEN (UNIT=19, STATUS="UNKNOWN" . FILE="talO . 300") 


C**  ** 
C**  CALL  SUBROUTINES  TO  SET  UP  ORIGINAL  POSITIONS  ** 
C**  AND  ORIGINAL  MOMENTUMS  ** 
C**  ** 


CALL  POS(X,Y,T) 
ORIGMOM  =  100. 

CALL  MOM (XMOM, YMOM, T) 
C 


oonnnnnnon  nnonnn  nn 


UNITS 


Mass  =  atomic  mass  units  (amu) 
Distance  =  Angstrom  (A) 

Time  =  picosecond  (psec) 


CALL  SUBROUTINE  TO  FIND  TOTAL  ENERGY  OF 
SYSTEM  AT  INITIAL  CONDITIONS 


CALL  TOTEN (X, Y, XMOM, YMOM, T) 

"DIVIDE"  IS  USED  TO  BREAK  THE  OUTPUT  INTO  "DIVIDE"  NUMBER  OF 
OUTPUT  FILES 

DIVIDE  =  10. 


DO  110  TIM  =  TINIT, TFINAL , TINC 
C 
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Q*  it  *  * 

C**  CALL  SUBROUTINES  TO  FIND  VELOCITIES  AND  MOMENTUM  DERIVATIVES  ** 
C**  ** 


IT  =  4 

CALL  WRITFILE  (X,  Y,  T,  IT,  TIM) 

CALL  LJ (X, Y,XMOM, YMOM,XVEL, YVEL,XMOMD, YMOMD,T) 


C**  USE  EULER  METHOD  TO  FIND  TEMPORARY  NEW  POSITIONS  AND  MOMENTUMS  ** 
C**  ** 


DO  120  I  =  1,N,1 

X2  (I)  =  X(I)  +  TINC  *  XVEL(I) 

Y2  (I)  =  Y(I)  +  TINC  *  YVEL(I) 

XMOM2  (I)  =  XMOM(I)  +  TINC  *  XMOMD(I) 
YMOM2  (I)  =  YMOM(I)  +  TINC  *YMOMD(I) 
120  CONTINUE 


C  ******************************************************************** 

C**  ** 

C**  CALL  SUBROUTINES  TO  FIND  VELOCITIES  AND  MOMENTUM  ** 

C**  DERIVATIVES  USING  ** 

C**  TEMPORARY  POSITIONS  AND  MOMENTUMS  FOR  USE  IN  MODIFIED  EULER  ** 
C **  METHOD  ** 

C**  ** 


CALL  LJ  (X2 , Y2 , XMOM2 , YMOM2 , XVEL2 , YVEL2 . XMOMD2 , YMOMD2 , T) 


C**  USE  MODIFIED  EULER  METHOD  TO  FIND  ACTUAL  NEW  POSITIONS  AND 

C**  METHODS 

C** 


non  on  on  o  ooooo  ooooo 


Y (I)  =  Y(I)  +  TINC2  *  (YVEL(I)  +  YVEL2(I)) 

XMOM(I)  =  XMOM(I)  +  TINC2  *  (XMOMD(I)  +  XM0MD2(I)) 

YMOM  (I)  =  YMOM  (I)  +  TINC2  *  (YMOMD(I)  +  YM0MD2  (I)  ) 

130  CONTINUE  I 


CHECK  TO  SEE  IE  TIM  =  ONE  OF  THE  DESIRED  OUTPUT  FLAGS 


DO  140  KONT  =  1,10,1 
CNT  =  FLOAT (KONT) 

IF  (ABS(TIM- (CNT*TFINAL/DIVIDE) ) -LT.1E-6)  THEN 
I  =  KONT  +9 

CALL  WRITFILE (X, Y, T, I , TIM) 

CALL  TOTEN (X, Y, XMOM, YMOM, T) 

CALL  STORE I LE (TIM, X, Y, XMOM, YMOM, T) 

END  IF 

140  CONTINUE 


MOVE  ATOM  IF  IT  HAS  GONE  OUTSIDE  OF  ORIGINAL  CRYSTAL  PARAMETERS 


DO  150  K  =  1 ,  N,  1 

IF  (Y  (K)  .  GT .  YMOVE)  Y  (K)  =  Y  (K)  -  YMOVE 

IF  (Y(K).LT.O.)  Y  (K)  =  Y  (K)  +  YMOVE 
IF  (X(K).LT.O.)  THEN 
XMOM  (K)  =  -XMOM  (K) 

X  (K)  =  -  X  (K) 

END  IF 

150  CONTINUE 
110  CONTINUE 

CALCULATE  THE  TOTAL  ENERGY  OF  THE  SYSTEM  AT  THE  END  OF  THE  PROGRAM 
CALL  TOTEN (X,Y, XMOM, YMOM, T) 

END 


SUBROUTINE  LJ  (X,  Y,  XMOM,  YMOM,  XVEL,  YVEL,  XMOMD,  YMOMD,  T) 


DIMENSION  X  ( * )  ,  Y  ( * )  ,  XMOM  ( * )  ,  YMOM  ( * ) 
DIMENSION  XVEL  (*)  ,YVEL(*)  ,  XMOMD  (*)  ,  YMOMD  (*) 
DIMENSION  T  (*) 

PARAMETER  (N  =  144) 

PARAMETER  (XORIGIN  =15.) 

PARAMETER  (SIZES  =  3.69) 

PARAMETER  (YMOVE  =  6.*SIZES) 

COORDINATE  DERIVATIVES 

DO  200  I  =  1 , N, 1 
C  REM  TO  REZERO  XMOMD  AND  YMOMD 
XMOMD  (I)  =  0. 

YMOMD (I)  =  0. 

C  FIND  THE  VELOCITY  FOR  EACH  ATOM 
IF (T (I) .LT.3.5)  THEN 


on  non  on 
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AMASS  =  1.008 
ELSE 

AMASS  =  35.45 
ENDIE 

YVEL(I)  =  YMOM(I)  /  AMASS 
200  XVEL  (I)  =  XMOM  (I)  /  AMASS 

C  ENLARGE  THE  CRYSTAL  TO  GET  STABILITY  FOR  THE  MONITERED  ATOMS 
DO  210  I  =  1,N,  1 
IF  (Y(I)  .GT.  0.5  *  YMOVE)  THEN 
Y(N  +  I)  =  Y(I)  -  YMOVE 
X  (N  +  I)  =  X(I) 

ELSE 

Y(N  +  I)  =  Y(I)  +  YMOVE 
X  (N  +  I)  =  X(I) 

END  IF 

210  CONTINUE 

C  I  HAVE  SET  UP  MOLECULES  ABOVE  AND  BELOW  THE  INITIAL  CRYSTAL 
C  THEN  GOING  THROUGH  THE  LIST  ONCE  CALCULATE  ALL  OF  THE  FORCES 
CALL  AL (X. Y, XMOM. YMOM, XVEL , YVEL , XMOMD , YMOMD , T) 

RETURN 

END 


SUBROUTINE  AL  (X.  Y,  XMOM,  YMOM,  XVEL ,  YVEL ,  XMOMD ,  YMOMD ,  T) 
DIMENSION  X(*)  ,Y  (’•'),  XMOM  (*),  YMOM  (*) 

DIMENSION  XVEL  (*)  ,YVEL(*)  ,  XMOMD  (*)  ,  YMOMD  (*) 
DIMENSION  T (*) 

PARAMETER  (N  =  144) 

PARAMETER  (SIZES  =  3.69) 

PARAMETER  (YMOVE  =  6.*SIZES) 

INTERMEDIATE  FUNCTIONS 


DO  300  I  =  l.N.l 
DO  310  J  =  1,  2*N  ,1 

FIND  THE  DISTANCE  BETWEEN  EACH  SET  OF  TWO  MOLECULES 
in  EACH  COORDINATE  DIRECTION  (X  &  Y) 

IF  (I.EQ.J)  GO  TO  310 
RX  =  X(I)  -  X(J) 

RY  =  Y(I)  -  Y  (J) 

C  FIND  TOTAL  DISTANCE  BETWEEN  EACH  SET  OF  TWO  MOLECULES 
R2  =  RX  *  RX  +  RY  *  RY 
R6  =  R2  *  R2  *  R2 
R  =  SQRT  (R2) 

C  FOR  LARGE  VALUE  OF  R  THE  FORCE  EXERTED  BY  MOLECULES  IS  SMALL 
C  THEREFORE  SAVE  TIME  BY  NOT  CALCULATING  THE  FORCE  IF  R  IS  TOO 
C  LARGE  (SKIP  BY  USING  IF  THEN  STATEMENT) 

IF (R.GT.l .E6)  GO  TO  310 

C  FIND  THE  LENNAED- JONES  FUNCTION  BETWEEN  EACH  SET  OF  TWO  MOLECULES 
C  REMEMBER  THAT  THE  MOLECULES  EXERT  EQUAL  AND  OPPOSITE  FORCES  ON  ANOTH 
IF (R.LT.O .0005*SIZES)  THEN 
WRITE(0, 320) 

320  FORMAT ("PARTICLES  ARE  GETTING  TOO  CLOSE  ERROR  HAS  OCCURRED") 
FPART  =  1000. 

ELSE 

IF  (  (T(I)  .LT.3.5)  .AND.  (T  (J)  .  LT .  3 . 5)  )  THEN 


onnn 


May  23  14:53  1986  deug.f  Page  7  A-7a 


C  A  amu  *  A* *8  /  psec**2  t 

A  =  15032.44  ; 

C  B  amu  *  A* *14  /  psec**2  ] 

B  =  1296.3691  I 

ELSE 

IF  (  (T(I)  .GT.3.5)  .AND.  (T(J)  .GT.  3.5)  )  THEN 
A  =  2998181.1 
B  =  92505858. 

ELSE 

A  =  368924.09 
B  =  788575.56 
END  IF 
ENDIF 

FPART  =6.  *  (2.  *  B  /  R6  -  A)  /  (R2*R6) 

ENDIF 

MOMENTUM  DERIVATIVES  * 

THE  RX  AND  RY  TERMS  ACCOUNT  FOR  THE  DIRECTION  OF  THE  FORCE 
XVDSUB  =  FPART  *  RX 
YVDSUB  =  FPART  *  RY 
C  SUM  ALL  THE  INDIVIDUAL  FORCES  FOR  EACH  MOLECULE 
XMOMD(I)  =  XMOMD(I)  +  XVDSUB 
YMOMD(I)  =  YMOMD(I)  +  YVDSUB 
310  CONTINUE 
300  CONTINUE 
RETURN 
END 

SUBROUTINE  MOM  (XMOM,  YMOM,  T) 

DIMENSION  XMOM (*)  ,  YMOM  (*) 

DIMENSION  T  (*) 

COMMON /B 1 / OR I GMOM 
PARAMETER  (N  =  144) 

C  SET  ORIGINAL  MOMENTUMS  TO  0. 

DO  400  I  =  1 ,  N,  1 
XMOM  (I)  =  0. 

YMOM (I)  =  0. 

400  CONTINUE 

C  GIVE  FIRST  FOUR  ATOMS  IN  EACH  ROW  INITIAL  MOMENTUM 
DO  410  I  =  1,6,1 
NSTART  =  24  *  (1-1)  +  1 
NEND  =  NSTART  +  3 
DO  420  J  =  NSTART, NEND, 1 

Y  =  RAND  (J) 

Y  =  (Y  -  0.5)  *  50. 

YMOM (J)  =  Y 

IF (T (J) . LT .3.5)  THEN 
XMOM (J)  =  ORIGMOM 
ELSE 

XMOM ( J)  =  35  *  ORIGMOM 
ENDIF 

420  CONTINUE 
410  CONTINUE 
RETURN 
END 

SUBROUTINE  TOTEN  (X,  Y,  XMOM,  YMOM,  T) 
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DIMENSION  T (*) 

DIMENSION  X(*)  ,Y(*) 

DIMENSION  XMOM(*)  ,YMOM(*) 
PARAMETER  (N  =  144) 

PARAMETER  (SIZES  =  3.69) 
PARAMETER  (YMOVE  =  6.  *SIZES) 

DO  500  I  =  1,N,  1 

IF  (Y(I) .GT.  0.5  *  YMOVE)  THEN 


I) 

I) 


500 


C  B 


C 

C 


Y  (N 
X  (N  + 
ELSE 
Y(N  + 

X(N  + 

END  IE 
CONTINUE 
ESYS  =  0, 
just  =  0 
DO  510  I 


=  Y(I) 

=  X(I) 


-  YMOVE 


I) 

I) 


=  Y(I) 
=  X(I) 


+  YMOVE 


YMOM(I)  *YMOM (I)  ) 


1,N,1 

XYMOMT2  =  (XMOM(I)  *XMOM(I) 

IF  (T  (I)  . LT .3.5)  THEN 
AMASS  =  1.008 
ELSE 

AMASS  =  35.45 
END  IF 

EKINTOT  =  XYMOMT2  /  2.  /  AMASS 
FTOT  =  0. 

DO  520  J  =  (1+1)  ,2*N,1 
RX  =  X(I)  -  X(J) 

RY  =  Y  (I)  -Y  (J) 

R2  =  RX*RX  +  RY*RY 

R6  =  R2  *  R2  *  R2 

IF  (R6.GT.1.E6)  GO  TO  520 

IF(J.GT.N)  THEN 

J1  =  J  -  N 

ELSE 

J1  =  J 

END  IF 

IF  (  (T(I)  .LT.3.5)  .AND.  (T(J1)  .LT.  3 .5) 


C  A  amu  *  A* *8  /  psec**2 


)  THEN 


A  =  15032.44 

amu  *  A* *14  /  psec**2 


B  =  1296.3691 
ELSE 

IF (  (T (I) .GT.3.5) 
A  =  2998181.1 
B  =  92505858. 
ELSE 

A  =  368924.09 
B  =  788575.56 
END  IF 
END  IF 


.AND.  (T(J1)  .GT.  3.5)  )  THEN 


CALL  STORFILE  (TIM,  X,  Y,  XMOM,  YMOM,  T) 
FPART  =  (B  /  (R6  *  R6) )  -  (A  /  R6) 
IF(J.GT.N)  FPART  =  0.5  *  FPART 


•*« A'". 
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FTOT  =  FTOT  +  FPART 
520  CONTINUE 

ESYS  =  ESYS  +  EKINTOT  +  FTOT 
510  CONTINUE 

WRITE (1, 530)  ESYS 

530  FORMAT (IX, "THE  SYSTEM'S  ENERGY  IS  ",E12.4) 
RETURN 
END 


SUBROUTINE  STORE ILE (TIM , X , Y , XMOM , YMOM , T) 
DIMENSION  X(*)  ,Y(*)  ,XMOM(*)  ,YMOM(*) 

DIMENSION  T  (*) 

PARAMETER  (N  =  144) 

REWIND  2 
WRITE (2 , 600) 

600  FORMAT (IX, "FILE  CONTAINING  TIM, I , T, X, Y, XMOM, 
WRITE (2, 610) TIM 
610  FORMAT ( IX, F7. 4) 

DO  620  I  =  1,N, 1 

WRITE  (2, 630)  I, T (I)  ,X(I)  ,  Y (I)  ,  XMOM (I)  ,YMOM(I) 
630  FORMAT (IX, 14, ",  ", F4 . 1, ” , F10 . 6 , , F10 . 6, " , 
620  CONTINUE 
RETURN 
END 

SUBROUTINE  POS(X.Y,T) 

DIMENSION  T  (*) 

DIMENSION  X(*)  ,  Y (*) 

C  set  up  array  for  type  of  molecule 
C* ************************************************ 


and  YMOM") 


", F13 .6, ",  ’ 


c** 

★  * 

c** 

KEY 

** 

c** 

*  * 

c** 

H-H  1 

★  * 

c** 

H  2 

** 

c** 

(H)  -Cl  3 

*  * 

c** 

H-  (Cl)  4 

*  * 

c** 

Cl  5 

★  * 

c** 

Cl -Cl  6 

*  * 

c** 

★  * 

c** 

"k  * 

PARAMETER  (N  =  144) 

PARAMETER (XORIGIN  =15.) 

PARAMETER (SIZES  =  3.69) 

PARAMETER (LENGTH  =  10.11411) 

C 

C 

C 

c** 

** 

c** 

Set  up  Original  Positions 

★  ★ 

c** 

*  * 

c** 

and 

★  ★ 

c** 

*  * 

c** 

Label  the  Molecules  by  Type 

*  * 
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Q*  ********************************************************  * 

c ** 

C**  DISTHH  and  DISTCLCL  ARE  INTRAMOLECULAR  DISTANCES 

C** 

C**  DIS1HCL  IS  AN  INTERMOLECULAR  DISTANCE 

C** 

C**  DISTCENT  IS  THE  DISTANCE  REQUIRED  TO  ALIGN 

C**  THE  CENTER  OE  THE  BONDS  BETWEEN  ROWS 

C** 

C**  All  distances  are  in  Angstroms 

C** 

Q* ********************************************************* 


********** 

£ 

*  ★ 

*  * 

*  * 

■ 

*  * 

•  *; 

*  * 

*  * 

„  • 

*  * 

£ 

*  * 

*  * 

*★ 

********** 

V. 

C 

C 

DISTHH  =  .74611 
DISTHCL  =3.69 
DISTCLCL  =  1.988 
DISTCENT  =  .620945 
K  =  -3 
J  =  0 

DO  700  YPOS  =  .5*SIZES, 6 . 2*SIZES, SIZES 

DO  710  XPOS  =  XORIGIN,  XORIGIN  +  5.1  *  LENGTH,  LENGTH 

K  =  K  +  4 

IE(J.EQ.O)  THEN 

X  (K)  =  XPOS 

X (K+2)  =  XPOS  +  DISTHH 
X (K+l)  =  XPOS  +  DISTHH  +  DISTHCL 
X (K+3)  =  XPOS  +  DISTHH  +  DISTHCL  +  DISTCLCL 
T  (K)  =1. 

T  (K+2)  =  1. 

T  (K+l)  =  6. 

T  (K+3)  =  6. 

ELSE 

X (K)  =  XPOS  -  DISTCENT 

X (K+2)  =  XPOS  -  DISTCENT  +  DISTCLCL 

X (K+l)  =  XPOS  -  DISTCENT  +  DISTCLCL  +  DISTHCL 

X (K+3)  =  XPOS  -  DISTCENT  +  DISTCLCL  +  DISTHCL  +  DISTHH 

T(K)  =6. 

T  (K+2)  =  6. 

T  (K+l)  =  1. 

T  (K+3)  =  1. 

END  IF 

Y  (K)  =  YPOS 

Y (K+2)  =  YPOS 
Y (K+l)  =  YPOS 

Y  (K+3)  =  YPOS 
710  CONTINUE 

IF(J.EQ.O)  THEN 

J  =  1 

ELSE 

J  =  0 

END  IF 

700  CONTINUE 

C  MOVE  FIRST  FOUR  ATOMS  IN  EACH  ROW  BACK 
DO  410  I  =  1,6,1 
NSTART  =  24  *  (1-1)  +  1 


.v 


* 
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NEND  =  NSTART  +  3 
DO  420  J  =  NSTART , NEND , 1 
X(J)  =  X(J)  -  5. 

420  CONTINUE 
410  CONTINUE 
1  =  3 

CALL  WRITFILE  (X.  Y, T, I ,  TIM) 

RETURN 

END 

SUBROUTINE  WRITFILE  (X,  Y,T,  I, TIM) 
DIMENSION  X(*)  ,Y(*)  ,T(*) 
PARAMETER  (N  =  144) 

REWIND  I 
WRITE (I, 800) 

800  FORMAT (IX, "display  world;") 
WRITE (I, 810) 
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WRITE  (I,  810) 

FORMAT ("  ") 

WRITE (I , 820) 

FORMAT ( IX , "  label : =begin_ 
WRITE (I , 830) 

FORMAT (IX,  ‘  charact 

WRITE  (I.  840) 

FORMAT ( IX , "  charact 

J  =  I  -  9 

WRITE  (1 , 850)  TIM 

FORMAT  (IX,"  charact 

WRITE (I, 860) 

FORMAT ( IX , "  encLstructure ; " ) 
WRITE (I, 810) 

WRITE (I, 810) 


label : =begin_structure") 

character  scale  5,5;") 

characters  25,50  'Modified  Euler  Method';" 


characters  25,40  'at  ”,F6.4,"  seconds’;") 


DO  870  J  =  l.N, 1 
IF  (T(J)  .LT.3.5) 
IF  (T(J)  .GT.3.5) 
RADIUS  =  RDIS  * 


RDIS  =  .37 
RDIS  =  .99 
0.7010678 


PY1  =  Y(J)  +  RDIS 
PX2  =  X(J)  -  RADIUS 
PY2  =  Y(J)  +  RADIUS 
PX3  =  X(J)  -  RDIS 
PX4  =  X(J)  -  RADIUS 
PY4  =  Y (J)  -  RADIUS 
PY5  =  Y(J)  -  RDIS 
PX6  =  X(J)  +  RADIUS 
PY6  =  Y(J)  "  RADIUS 
PX7  =  X(J)  +  RDIS 
PX8  =  X(J)  +  RADIUS 
PY8  =  Y(J)  +  RADIUS 
K1  =  J  +  100 
WRITE  (I, 880)  K1 

FORMAT (IX, "c", 13, " : =vector_list  n  =  20") 
WRITE  (I,  890)  X(J)  ,  PY1 
WRITE  (1 , 890)  PX2 ,  PY2 
WRITE  (1 , 890)  PX3,  Y  (J) 

WRITE  (1,890)  PX4.PY4 


P* 
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WRITE  (1 , 890)  X  (J)  , PY5 
WRITE  (1 , 890)  PX6 ,  PY6 
WRITE  (1 , 890) PX7, Y  (J) 

WRITE (I , 890) PX8 , PY8 
WRITE  (1 , 890)  X  (J)  .  PY1 
FORMAT ("  " , F9 .4, "  " 
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.  F9 . 4) 


WRITE  (I , 900) 

FORMAT  ("  ;") 

WRITE  (1 , 910) 

WRITE  (I,  910) 

FORMAT  ("  ") 

CONTINUE 
WRITE (1 , 920) 

FORMAT  (IX, "world: =begin_structure") 

WRITE  (1 ,930) 

FORMAT (IX, "  window  x=0:150  y=-20:130*") 

DO  940  K  =  1,  N,  1 
K1  =  K  +  100 
IF  (T (K)  .EQ.l.)  J  =  120 
IF  (T  (K)  .  EQ .  2  . )  J  =  168 
IF  (T  (K)  .  EQ .  3 . )  J  =  359 
IF  (T  (K)  .  EQ .  4 . )  J  =  359 
IF  (T  (K)  .  EQ .  5 . )  J  =  216 
IF  (T(K)  .  EQ .  6 . )  J  =  264 
WRITE  (1 , 950)  J,  K1 
FORMAT (IX. "  set  color  ",I3,",1  applied  to  c",I3, 

CONTINUE 
WRITE (I, 955) 

FORMAT (IX,"  instance  of  label;") 

WRITE (I, 960) 

FORMAT (IX, "  endLstructure; ") 

RETURN 
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SUBROUTINE  RECALE  (X,  Y,  XMOM,  YMOM,  T,  TIMA) 
DIMENSION  X  (*)  ,  Y(*)  ,  XMQM(*)  ,  YMQM(*) 

DIMENSION  T  (*) 

PARAMETER  (N  =  144) 

READ (8, 300) TIMA 
300  FORMAT  (F8. 5) 

DO  200  I  =  1,N,  1 

READ  (8, 400)  A,  INT,  T  (I)  ,C,X(I)  ,  B,  Y  (I)  ,  D,  XMOM  (I)  ,E 
400  FORMAT (F3.1, 14, 2F4. 1, F10 .6, F3 . 1, F10 .6. F3 . 1, F13 . 
200  CONTINUE 

WRITE  (0,  3)  TIMA 
3  FORMAT ( lx , f 8. 5) 

RETURN 


,  YMOM  (I) 

6, F3 . 1, F13 . 6) 
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label 0 : =begin_structure 

character  scale  3,3; 

characters  11,50  'Modified  Euler  Method'; 
characters  11,40  'at  0-0000  picoseconds'; 
end_s tructur e ; 


cOlOl :=vector_list  n  =  20 
5.0000,  2.2150 

4.7406,  2.1044 

4.6300,  1.8450 

4.7406,  1.5856 

5.0000,  1.4750 

5.2594,  1.5856 

5.3700,  1.8450 

5.2594,  2.1044 

5.0000,  2.2150 


C0102 :=vector_list  n  =  20 

9.4361,  2.8350 

8.7421,  2.5391 

8.4461,  1.8450 

8.7421,  1.1509 

9.4361,  0.8550 

10.1302,  1.1509 

10.4261,  1.8450 

10.1302,  2.5391 

9.4361,  2.8350 


c0103:=vector_list  n  =  20 

5.7461,  2.2150 

5.4867,  2.1044 

5.3761,  1.8450 

5.4867,  1.5856 

5.7461,  1.4750 

6.0055,  1.5856 

6.1161,  1.8450 

6.0055,  2.1044 

5.7461,  2.2150 


c0104:=vector_list  n  =  20 
11.4241,  2.8350 

10.7301,  2.5391 

10.4341,  1.8450 

10.7301,  1.1509 

11.4241,  0.8550 

12.1182,  1.1509 

12.4141,  1.8450 
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12 . 1182, 

11.4241, 


2.5391 

2.8350 


c0105 :=vector_list  n  =  20 
25.0000,  2.2150 

24.7406,  2.1044 

24.6300,  1.8450 

24.7406,  1.5856 

25.0000,  1.4750 

25.2594,  1.5856 

25.3700,  1.8450 

25.2594,  2.1044 

25.0000,  2.2150 


c0106 :=vector_list  n  =  20 

29.4361,  2.8350 

28.7421,  2.5391 

28.4461,  1.8450 

28.7421,  1.1509 

29.4361,  0.8550 

30.1302,  1.1509 

30.4261,  1.8450 

30.1302,  2.5391 

29.4361,  2.8350 


c0107 : =vector_l 1st  n  =  20 

25.7461,  2.2150 

25.4867,  2.1044 

25.3761,  1.8450 

25.4867,  1.5856 

25.7461,  1.4750 

26.0055,  1.5856 

26.1161,  1.8450 

26.0055,  2.1044 

25.7461,  2.2150 


c0108 :=vector_list  n  =  20 


31.4241, 

30.7301, 

30.4341, 

30.7301, 

31.4241, 

32.1182, 
32.4141, 

32.1182, 

31.4241. 


2.8350 

2.5391 

1.8450 

1.1509 

0.8550 

1.1509 

1.8450 

2.5391 

2.8350 
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c0109:=vector_list  n  =  20 
35 . 0000 ,  2.2150 

34.7406,  2.1044 

34.6300,  1.8450 

34.7406,  1.5856 

35.0000,  1.4750 

35.2594,  1.5856 

35.3700,  1.8450 

35.2594,  2.1044 

35.0000,  2.2150 


cOHO :=vector_list  n  =  20 


39.4361, 

38.7421, 

38.4461, 

38.7421, 

39.4361, 

40.1302, 

40.4261, 

40.1302, 

39.4361, 


2.8350 
2.5391 
1 . 8450 
1.1509 
0.8550 
1.1509 
1.8450 
2.5391 
2.8350 


cOlll :=vector_list  n  =  20 

35.7461,  2.2150 

35.4867,  2.1044 

35.3761,  1.8450 

35.4867,  1.5856 

35.7461,  1.4750 

36.0055,  1.5856 

36.1161,  1.8450 

36.0055,  2.1044 

35.7461,  2.2150 


c0112 :=vector_list  n  =  20 


41.4241, 

40.7300, 

40.4341, 

40.7300, 

41.4241, 

42.1182, 

42.4141, 

42.1182, 

41.4241, 


2.8350 

2.5391 

1.8450 

1.1509 

0.8550 

1.1509 

1.8450 

2.5391 

2.8350 


VAV-.V 


Y/A  VSl  3X* 
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c0244:=vector_list  n  =  20 

70 . 8032 , 

20.6650 

70.5438, 

20.5544 

70.4332, 

20.2950 

70.5438, 

20.0356 

70.8032, 

19.9250 

71.0626. 

20.0356 

71.1732. 

20.2950 

71.0626, 

20.5544 

70.8032. 

20.6650 

frameO : =begin_structure 


window  x=-30 : 120 

y=-60:90; 

set 

color 

120,1 

applied 

to 

C0101 

set 

color 

264,1 

applied 

to 

C0102 

set 

color 

120,1 

applied 

to 

C0103 

set 

color 

264,1 

applied 

to 

C0104 

set 

color 

120,1 

applied 

to 

C0105 

set 

color 

264.1 

applied 

to 

C0106 

set 

color 

120,1 

applied 

to 

C0107 

set 

color 

264,1 

applied 

to 

C0108 

set 

color 

120,1 

applied 

to 

C0109 

set 

color 

264,1 

applied 

to 

cOHO 

set 

color 

120,1 

applied 

to 

com 

set 

color 

264,1 

applied 

to 

C0112 

set 

color 

120,1 

applied 

to 

C0113 

set 

color 

264,1 

applied 

to 

C0114 

set 

color 

120,1 

applied 

to 

C0115 

set 

color 

264,1 

applied 

to 

C0116 

set 

color 

120,1 

applied 

to 

C0117 

set 

color 

264,1 

applied 

to 

C0118 

set 

color 

120,1 

applied 

to 

C0119 

set 

color 

264,1 

applied 

to 

C0120 

set 

color 

120,1 

applied 

to 

C0121 

set 

color 

264,1 

applied 

to 

C0122 

set 

color 

120,1 

applied 

to 

C0123 

set 

color 

264,1 

applied 

to 

C0124 
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set  color  264, 
set  color  120, 
set  color  264, 
set  color  120, 
set  color  264, 
set  color  120, 
set  color  264, 
set  color  120, 
set  color  264, 
set  color  120, 
set  color  264, 
set  color  120, 
set  color  264, 
set  color  120, 
set  color  264, 
set  color  120, 
set  color  264, 
set  color  120, 
set  color  264, 
set  color  120, 
set  color  264, 
set  color  120, 
set  color  264, 
set  color  120, 
set  color  120, 
set  color  264, 
set  color  120, 
set  color  264, 
set  color  120, 
set  color  264, 
set  color  120, 
set  color  264, 
set  color  120, 
set  color  264, 
set  color  120, 
set  color  264, 
set  color  120, 
set  color  264, 
set  color  120, 
set  color  264, 
set  color  120, 
set  color  264, 
set  color  120, 
set  color  264, 
set  color  120, 
set  color  264, 
set  color  120, 
set  color  264, 
set  color  264, 
set  color  120, 
set  color  264, 
set  color  120, 
set  color  264, 
set  color  120, 
set  color  264, 
set  color  120, 
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1  applied  to  c0125; 
1  applied  to  c0126; 
1  applied  to  c0127; 
1  applied  to  c0128; 
1  applied  to  c0129; 
1  applied  to  c0130; 
1  applied  to  c0131; 
1  applied  to  c0132; 
1  applied  to  c0133; 
1  applied  to  c0134; 
1  applied  to  c0135; 
1  applied  to  c0136; 
1  applied  to  c0137; 
1  applied  to  c0138; 
1  applied  to  c0139; 
1  applied  to  c0140; 
1  applied  to  c0141; 
1  applied  to  c0142; 
1  applied  to  c0143; 
1  applied  to  c0144; 
1  applied  to  c0145; 
1  applied  to  c0146; 
1  applied  to  c0147; 
1  applied  to  c0148; 
1  applied  to  c0149; 
1  applied  to  c0150; 
1  applied  to  c0151; 
1  applied  to  c0152; 
1  applied  to  c0153; 
1  applied  to  c0154; 
1  applied  to  c0155; 
1  applied  to  c0156; 
1  applied  to  c0157; 
1  applied  to  c0158; 
1  applied  to  c0159; 
1  applied  to  c0160; 
1  applied  to  c0161; 
1  applied  to  c0162; 
1  applied  to  c0163; 
1  applied  to  c0164; 
1  applied  to  c0165; 
1  applied  to  c0166; 
1  applied  to  c0167; 
1  applied  to  c0168; 
1  applied  to  c0169; 
1  applied  to  c0170; 
1  applied  to  c0171; 
1  applied  to  c0172; 
1  applied  to  c0173; 
1  applied  to  c0174; 
1, applied  to  c0175; 
1  applied  to  c0176; 
1  applied  to  c0177; 
1  applied  to  c0178; 
1  applied  to  c0179; 
1  applied  to  c0180; 
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set 

color 

264, 1 

set 

color 

120,1 

set 

color 

264,1 

set 

color 

120,1 

set 

color 

264, 1 

set 

color 

120,1 

set 

color 

264,1 

set 

color 

120,1 

set 

color 

264,1 

set 

color 

120,1 

set 

color 

264,1 

set 

color 

120,1 

set 

color 

264,1 

set 

color 

120,1 

set 

color 

264,1 

set 

color 

120,1 

set 

color 

120,1 

set 

color 

264,1 

set 

color 

120,1 

set 

color 

264.1 

set 

color 

120,1 

set 

color 

264,1 

set 

color 

120.1 

set 

color 

264,1 

set 

color 

120,1 

set 

color 

264,1 

set 

color 

120,1 

set 

color 

264,1 

set 

color 

120,1 

set 

color 

264, 1 

set 

color 

120,1 

set 

color 

264,1 

set 

color 

120,1 

set 

color 

264, 1 

set 

color 

120,1 

set 

color 

264,1 

set 

color 

120,1 

set 

color 

264,1 

set 

color 

120,1 

set 

color 

264, 1 

set 

color 

264, 1 

set 

color 

120,1 

set 

color 

264, 1 

set 

color 

120,1 

set 

color 

264, 1 

set 

color 

120,1 

set 

color 

264,1 

set 

color 

120,1 

set 

color 

264,  1 

set 

color 

120, 1 

set 

color 

264,1 

set 

color 

120,1 

set 

color 

264,1 

set 

color 

120,1 

set 

color 

264,1 

set 

color 

120,1 

applied  to  c0181; 
applied  to  c0182 
applied  to  c0183 
applied  to  c0184 
applied  to  c0185 
applied  to  c0186 
applied  to  c0187 
applied  to  c0188 
applied  to  c0189 
applied  to  c0190 
applied  to  c0191 
applied  to  c0192 
applied  to  c0193 
applied  to  c0194 
applied  to  c0195 
applied  to  c0196 
applied  to  c0197 
applied  to  c0198 
applied  to  c0199 
applied  to  c0200 
applied  to  c0201 
applied  to  c0202 
applied  to  c0203 
applied  to  c0204 
applied  to  c0205 
applied  to  c0206 
applied  to  c0207 
applied  to  c0208 
applied  to  c0209 
applied  to  c0210 
applied  to  c0211 
applied  to  c0212 
applied  to  c0213 
applied  to  c0214 
applied  to  c0215 
allied  to  c0216 
applied  to  c0217 
applied  to  c0218 
applied  to  c0219 
applied  to  c0220 
applied  to  c0221 
applied  to  c0222 
applied  to  c0223 
applied  to  c0  224 
applied  to  c0225 
applied  to  c0226 
applied  to  c0227 
applied  to  c0228 
applied  to  c0229 
applied  to  c0230 
applied  to  c0231 
applied  to  c0232 
applied  to  c0233 
applied  to  c0234 
applied  to  c0235 
applied  to  c0236 
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set  color  264,1  applied  to  c0237; 
set  color  120,1  applied  to  c0238; 
set  color  264,1  applied  to  c0239; 
set  color  120,1  applied  to  c0240; 
set  color  264,1  applied  to  c0241; 
set  color  120,1  applied  to  c0242; 
set  color  264,1  applied  to  c0243; 
set  color  120,1  applied  to  cQ244; 
instance  of  labelO; 
end_structur e ; 


